package gbench.sandbox.matlib;

import org.junit.jupiter.api.Test;

import java.text.NumberFormat;
import gbench.common.matlib.MatlibCanvas.NVec;
import gbench.common.matlib.MatlibCanvas.Seq;
import gbench.common.matlib.MatlibCanvas.XVec;

import static java.lang.Math.random;

import static gbench.common.matlib.MatlibCanvas.println;
import static gbench.common.matlib.data.Pipe.DUMMY;

public class JunitRegression {
    
    /**
     * 多元回归分析
     */
    @Test
    public void regression() {
        final var nf = NumberFormat.getInstance(); // 数字格式化函数
        
        final var volume = 300; // 数据数量
        final var num_variables = 10; // 变量数量,需要注意 一定要整除volume
        final var mx_pipe = DUMMY(Double[][].class).bind(XVec::OF).bind(xv -> xv.fmap(NVec::new)); // 源数据的加工管道, 管道类型为:Pipe<Double[][], XVec<NVec>>
        final var mx = NVec.of(volume, i -> random()).reshape2(volume/num_variables, num_variables); // 把一条一维数据向量reshape2转换成二维数组形式的数据结构
        
        final var x_train = mx_pipe.evaluate(mx); // 转转 换成 XVec<NVec> 类型的数据结构
        println("\nx_train 训练数据x矩阵:\n", x_train.toString(e -> e + "\n"));
        
        final var x_train_t = XVec.OF(Seq.transpose(mx)).fmap(NVec::new); // x_train的转置 也是 XVec<NVec> 类型的结构
        println("\nx_train_t 训练数据x矩阵的转置 :\n", x_train_t.toString(e -> e + "\n")); // 转置是为了方便取列数据 以方便计算：error.dot(x_train_t.get(j)) 即j索引方向的梯度分量
        
        final var actual_thetas = NVec.of(num_variables, i -> i + 1);// 估计的参数
        println("actual_thetas 训练数据参数实际值:\n", actual_thetas.toString(e -> e + "\n"));
        
        final var y_train = x_train.fmap(actual_thetas::dot).mutate(NVec::new); // 标准结果
        println("\ny_train 训练数据实际值:\n", y_train.toString(e -> e + "\n"));

        final var thetas = NVec.of(num_variables, i -> 0); // 估计的参数
        println("参数thetas的估计值（初始状态,默认为0):" + thetas.toString(nf::format));
        
        final var error_pipe = DUMMY((XVec<NVec>) null) // 输入一条观测行向量组：XVec<NVec> 即源数据矩阵 x_train , 管道类型为: Pipe<XVec<NVec>, NVec>
            .bind(xvec -> xvec.fmap(thetas::dot).mutate(NVec::new)) // x_i_j*theta_j -> y_i
            .bind(ynvec -> ynvec.sub(y_train)); // 残差计算管道 ynvec:[y_i]
        
        final var rate = 0.001; // 学习率-步长
        final var eps = 1e-10d; // 误差容许范围：最大值的误差累计和
        final var error = NVec.of(x_train.size(), i -> 0); // 残差向量，可以采用error_pipe.evaluate(x_train);来模拟计算一下error向量长度，但该结果在还会循环体重新计算一次。故此处使用0向量做占位。
        
        var i = 0; // 循环变量
        var loss = Double.MAX_VALUE; // loss 是一个error向量函数，把error向量生成一个loss指标。以便做容忍限度eps的检测，即调整thetas使得loss达到容许范围
        // 计算残差 以及 相应的 损失函数 直至 达到 低于 容许范围 eps要求
        while ( (loss = error.assign(error_pipe.evaluate(x_train)).fmapN(e -> e * e).mean()) >= eps ) { // LOSS = MSE（Mean Square Error ）均方误差
            println(i++, "#", loss); // 显示残差大小
            // 注意 - rate 是由 error_pipe 的 ynvec-y_train来定义的。符号(-)代表与 ynvec的增加方向（梯度方向)相反,即 负梯度方向
            thetas.update((j, theta_j) -> theta_j - rate * 2 * error.dot(x_train_t.get(j))); // thetas 系数重新计算
        } // while
        
        println("参数thetas的估计值（回归结果）:" + thetas.toString(nf::format));
        println("参数thetas估计误差 (thetas-actual_thetas):\n",thetas.sub(actual_thetas));
    }

}
